

#copyright 2009 Kevin Daley
##This file is a part of Nyst.
#Nyst is free software: you can redistribute it and/or modify
 #   it under the terms of the GNU Lesser General Public License as published by
  #  the Free Software Foundation, either version 3 of the License, or
   # (at your option) any later version.

   # This program is distributed in the hope that it will be useful,
   # but WITHOUT ANY WARRANTY; without even the implied warranty of
   # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   # GNU Lesser General Public License for more details.

   # You should have received a copy of the GNU Lesser General Public License
   # along with this program.  If not, see <http://www.gnu.org/licenses/>.

#Kevin Daley, Lead Developer.  still.horse@gmail.com.

from cmath import *
import numpy
from numpy import *
from random import * 
from scipy.ndimage import *
from scipy import *
from quadriterate import *
class network{
   def __init__(self, resolution, spectrum, number, cell_width, mass)
	self.network=numpy.sum(number*spectrum,axis=0)
	self.res_x=resolution[0]
	self.res_y=resolution[1]
	self.spectrum=spectrum
	self.number=number
	self.cellwidth=cell_width
	self.lambda2=[numpy.sum(number[k])/(res_x*res_y*d**2) for k in spectrum]
	self.mass=mass
	self.lambda1=rate_param[0]
	self.mesh_refinement_scale=resolution[2]
	a0=self.network
	Nx=self.res_x
	Ny=self.res_y
	spectrum=self.spectrum
	nk=self.number
	d=self.cellwidth
	lamb2a=self.lambda2
	lamba=self.lambda1
	finer=self.mesh_refinement_scale
        self.counter=0
        self.J=zeros((len(spectrum),Nx, Ny))
        t=0
	for m in xrange(len(spectrum)):
                        self.J[m]=array(Julia(W=Nx, H=Ny, c=a0[m]).image.getdata()).reshape(Nx,Ny)
   def evolve(self, loop, dt):
	self.timestep=dt
        a0=self.network
	Nx=self.res_x
	Ny=self.res_y
	spectrum=self.spectrum
	nk=self.number
	d=self.cellwidth
	dt=self.timestep
	lamb2a=self.lambda2
	lamba=self.lambda1
	finer=self.mesh_refinement_scale
        self.counter=0
        self.J=zeros((len(spectrum),Nx, Ny))
        t=0
        
	if(1):
                        for indek in xrange(len(spectrum)):
                       		k=indek
                                n=nk[indek]
                                julia=J[indek]
                                for i in xrange(0,Nx):
                                        for j in xrange(0,Ny):
                                                if julia[i,j]>0 or (i % finer==0 and j % finer==0):
                                                          for g in xrange(10):
                                                                  lamb2=lamb2a[indek]
                                                                  lamb=lamba[indek]
                                                                  n[i,j]-=int(abs((0.0+1.0j)*(i*d*lamb2/(1*k)+i*d*(1*k)*pow(lamb2,2)-(i*d*lamb/(k)-i*d*(k)*pow(lamb,2)))))
                                                                  network[i,j]-=self.mass*exp((0.0+1.0j)*dt*k)*(0.0+1.0j)*(i*d*lamb2/(1*k)+i*d*(1*k)*pow(lamb2,2)-(i*d/(lamb*k)-i*d*(k)*pow(lamb,2)))*k
                                                                  lamb2a[indek]-=dt*(0.0+1.0j)*(i*d*lamb2/(1*k)+i*d*(1*k)*pow(lamb2,2))/(Nx*Ny*d**2)
                                                                  lamba[indek]-=dt*(0.0+1.0j)*(i*d*lamb/(k)-i*d*(k)*pow(lamb,2))/(Nx*Ny*d**2)
                                t+=dt    
                                nmax=numpy.max(n)         
                                if(self.counter % 30):                               
                                        
                                        grad_levelset=gaussian_filter(n,nmax)
                                        for i in xrange(0,Nx):
                                                for j in xrange(1, Ny):
                                                        if(j+1 < Ny and j-1 >= 0):
                                                                grad_levelset[i,j]=((grad_levelset[i,(j+1)]-grad_levelset[i,j])*d/4+(grad_levelset[i,(j-1)]-grad_levelset[i,(j+1)])*d/4)
                                        for i in xrange(0,Nx):
                                                for j in xrange(0,Ny):
                                                       if julia[i,j]==0 and (i % finer !=0 or j % finer !=0):
                                                               for g in xrange(6):
                                                                        n[i,j]-=int(abs((0.0+1.0j)*grad_levelset[i,j]))
                                                                        network[i,j]-=self.mass*exp((0.0+1.0j)*dt*k)*(0.0+1.0j)*grad_levelset[i,j]*k
                                                                        lamb2a[indek]-=(0.0+1.0j)*dt*grad_levelset[i,j]/(Nx*Ny*d**2)
                                                                        lamba[indek]-=(0.0+1.0j)*dt*grad_levelset[i,j]/(Nx*Ny*d**2)
                                                                        nk[indek]=n
				self.counter+=1
	
					
       
                





